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SEMI-LAGRANGIAN DISCONTINUOUS GALERKIN SCHEMES 
FOR SOME FIRST- AND SECOND-ORDER PARTIAL 
DIFFERENTIAL EQUATIONS 

OLIVIER BOKANOWSKI AND GIOREVINUS SIMARMATA 


Abstract. Explicit, unconditionally stable, high-order schemes for the ap¬ 
proximation of some first- and second-order linear, time-dependent partial dif¬ 
ferential equations (PDEs) are proposed. The schemes are based on a weak 
formulation of a semi-Lagrangian scheme using discontinuous Galerkin (DG) 
elements. It follows the ideas of the recent works of Grouseilles, Mehren- 
berger and Vecil (2010), Rossmanith and Seal (2011), for first-order equations, 
based on exact integration, quadrature rules, and splitting techniques for the 
treatment of two-dimensional PDEs. For second-order PDEs the idea of the 
scheme is a blending between weak Taylor approximations and projection on 
a DG basis. New and sharp error estimates are obtained for the fully dis¬ 
crete schemes and for variable coefficients. In particular we obtain high-order 
schemes, unconditionally stable and convergent, in the case of linear first-order 
PDEs, or linear second-order PDEs with constant coefficients. In the case of 
non-constant coefficients, we construct, in some particular cases, ’’almost” un¬ 
conditionally stable second-order schemes and give precise convergence results. 
The schemes are tested on several academic examples. 


1. Introduction 

In this paper we consider equations of the form 

Ut — -Trlaa'^D^u) + b ■ Vu + ru = 0, x € fi, t € (0,T), (I) 

where C is a box (with some boundary conditions on d^l), a (matrix), h 
(vector) and r (scalar) may be cc-dependent, at least Lipschitz continuous, together 
with an initial condition 

u(0, x) = Uo(x), X G fi, (2) 

with uq G The matrix a may be zero or positive semidefinite. Unless 

otherwise stated, we will in general assume periodic boundary conditions for 0 in 
order to avoid difficulties on the boundary. We will assume sufficient regularity on 
the data in order to have existence and uniqueness of weak solutions of 0-0 , and 
so that t -G u{t ,.) is in (^^([O, T], L^(f2)). 

We study and propose new semi-Lagrangian Discontinuous Galerkin schemes, 
also abbreviated ”SLDG” in this work, in order to approximate the solutions of 

0 - 0 . 

The semi-Lagrangian (SL) approach (see [13], or the textbook dH), is based 
on the approximation of the ’’method of characteristics”. By considering a weak 
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formulation of this principle, an explicit SLDG scheme is obtained. In the case 
of first-order PDEs with constant coefficient, our approach is based on a similar 
method as in the recent works of Crouseilles, Mehrenberger and Vecil [5] (for the 
Vlasov equation in plasma physics), Rossmanith and Seal |34j . However our ap¬ 
proach seems not to have been considered for variable coefficients. It is slightly 
different from the work of Qiu and Shu [3T] (see also Restelli et al [32]), where 
first a weak formulation of the PDE is considered, and then quadrature formulae 
are used (see also |33| for the original approach). Here we will furthermore intro¬ 
duce new SLDG schemes for second-order PDEs for which we prove stability and 
convergence results, and obtain higher-orders of accuracy when possible. 

First, in Section 2, we revisit the one-dimensional first-order advection equation 
with non-constant advection term b{x) (case cr = 0 in 0). We give a new un¬ 
conditional stability result, and convergence proof, extending similar results of [5|, 
[34] (or m) that was obtained for the case of a constant advection term. The 
unconditional stability property can be interesting when compared to a standard 
DG approach where a restrictive GEL condition must in general be considered [7]. 

Based on the operator construction for first-order advection, we then introduce, 
in Section 3, new schemes for linear second-order PDEs of type Q, in the form 
of explicit high-order SLDG schemes. These schemes are based, for the temporal 
discretization, on the use of ’’weak Taylor approximations”, see in particular the 
review book by Kloeden and Platen |2Q| (see also Kushner [21] and the review 
book by Kushner and Dupuis [22], Platen [30], Milstein [25], Talay [37], Pardoux 
and Talay [55], Menaldi [^, Camilli and Falcone [3], Milstein and Tretyakov [55] . 
m)- Such approximations where used by R. Ferretti in m as well as in Debrabant 
and Jakobsen m in the context of semi-Lagrangian schemes, using interpolation 
methods for the space variable. The problem of coupling such approximations with 
a spatial grid approximation, in particular using a high-order interpolation method, 
can be the stability and the convergence proof of the method. The Pi interpolation 
is known to be L°° stable, but it is only second-order accurate in space (for regular 
data). Some higher-order SL approximations have been proved to be stable (and 
convergent) for specific equations and under large GEL numbers (see [HllS]), or for 
some advection equations when the SL scheme can be reinterpreted in a weak form 
(we refer in particular to Ferretti’s work [ffldi])- 

The schemes of the present paper can be seen as projections of these approxima¬ 
tion on a discontinous Galerkin basis. We will in particular propose a second-order 
approximation (in time) corresponding to a Platen’s scheme |20[ Ghapter 14], but 
higher-order approximations (in time) could be obtained in the same way. The 
scheme will be proved to be also high-order in space, stable and convergent under 
a weak GEL condition (of the form Ax'^ < XAt for some constant A, where At and 
Ax denote the time and mesh steps). 

For the more simple case of second-order PDEs with constant coefficients, we 
also propose explicit and unconditionally stable schemes, high-order in space and 
up to third-order in time (higher-order can be obtained [2]). 

In section 4 we consider extensions to some linear two-dimensional PDEs. For 
first-order PDEs, we show how to combine the scheme with higher-order splitting 
techniques, like Strang’s splitting, but also Ruth’s third-order splitting [35], Forest’s 
fourth-order splitting m and Yoshida’s sixth-order splitting m (see also m 
and [41]). A splitting strategy to treat general second-order PDEs with constant 
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coefRcients is explained. The case of second-order PDEs with variable diffusion 


as 


coefficients is discussed but only treated in some specific cases (see Remark 4.4 
well as Example 7 and Example 8 of Sectio n . The general case will be treated 
in a forthcoming work (see however Remark 4.3). 

Einally in Section 5 we show the relevance of our approach on several academic 
numerical examples in one and two dimensions (using Cartesian meshes), including 
also a Black and Scholes PDE in mathematical finance. 

The advantage of the proposed schemes is that they combine the DG framework 
which allows high-order spatial accuracy and the potential of degree adaptivity, 
together with unconditional stability properties in the norm from the weak 
formulation of the semi-Lagrangian scheme. 

Note that our general strategy is to use a Cartesian grid, a particular one¬ 
dimensional advection scheme, and splitting techniques (for more standard Discon¬ 
tinuous Galerkin approaches, see for instance [5] or [H]). 

Ongoing works using the current approach concern the construction of higher- 
order schemes for general second-order PDEs [2], extensions to nonlinear PDEs 
arising from deterministic control [3] or from stochastic control. 

Acknowledgments. This work was partially supported by the EU under the 
7th Framework Programme Marie Curie Initial Training Network “FP7-PEOPLE- 
2010-ITN”, SADCO project, GA number 264735-SADCO. The first author also 
wishes to thank K. Debrabant for pointing out Platen’s works as well as an anony¬ 
mous referee for related references, which helped to simplify the presentation. We 
also thank D. Seal for useful comments and references. We are grateful to C.-W. 
Shu for pointing out problems in the preliminary version of the present work. 


2. Advection equation 


We first consider the semi-Lagrangian Discontinuous Galerkin scheme (SLDG 
for short) for the following one-dimensional first-order PDE, as in [9] 


{t,x) e (0,T) X D 
X G fl 


(3) 


\vt-\- b{x)vx = 0, 

[u(0,a;) = voix), 

where D = (x^nm, Xmax)^ together with periodic boundary conditions on D. 

In order to simplify the presentation and the proofs, we will assume that D = 
(0,1) and that 6 is a 1-periodic function. 

Let y = Px denote the solution of the differential equation 

y{t) = h{y(t)), tGR , , 

y{0) =x. 

We will also assume that b{-) is Lipschitz continuous. 

Let TV € N, iV > 1, At = ^ a time step and = nAt a time discretization. Let 

v^{x) := v{tn,x). 

By the method of characteristics, the solution of (|^ satisfies 


v-+^{x)=v-iyx{-At)). 


(5) 


Then we aim to obtain a fully discrete scheme. 

Let us consider a space discretization that is considered uniform for the sake of 
simplicity of presentation. Let Ax = for some integer M > 1, := 

Xmin -\-iAx, Vi = 0,..,M, and C := (aij_i ,a;j_|_i). Let fc € N. We define 14 as the 
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space of discontinuous-Galerkin elements on Q with polynomials of degree A:, that 
is: 


Vk = {veL\n,R):v\uePk,yi = 0,..,M-l} (6) 

where Pk denotes the set of polynomials of degree at most k. 

Remark 2.1. In the classical semi-Lagrangian approach, looking for vPfx), an 
approximation of vitn,x), a first ’’direct” iterative scheme for ([^ would be 

= K](y,.(-At)) (7) 

where [u"](a;) denotes some interpolation of the function u" at point x. We could 
take for instance a set of k + 1 values {x],f)a=[j,...,k in each interval 1^, and define 
the new polynomial such that := [u’^]{xl„ — bAt) for all a = 0,... ,k. 

However, given the discontinuities between the intervals li, this may lead to insta¬ 
bilities in the scheme m)- For instance, taking x]^ to be the Gauss quadrature 
points on each interval li is in general unstable (see Appendix^^ see also m)- 

Here we consider a Lagrange-Galerkin approach by taking the weak form of (§: 
for n = 0,..., iV — 1, find € 14 such that 

J u'’^'^^{x)(p{x)dx = j u’"'{yx(—^t))ip{x)dx, V(/3 S 14, (8) 

n n 

and for n = 0, find G Vk such that: 

J u^{x)ip{x)dx = J Vo{x)<p{x)dx, G 14 - (9) 

n Q 

From now on, we rewrite (|^ in the following abstract form : 

=TbAt{un. 

In the case of a constant coefficient b, yx{—At) = x — bAt, and u"(a: — bAt) is 
a piecewise constant polynomial. The integral Jj u"{x — bAt)ip{x)dx will have in 
general two regular parts. Each part involves a polynomial of degree at most 2k 
and the Gaussian quadrature rule with fc + 1 points is applied and is exact. At this 
stage the method is the same as in [^, or [^. Hence the new function can 
be computed by solving exactly Q. 

However, if b{x) is not a constant, x —>■ At)) is no more a piecewise 

polynomial. Therefore the computing procedure for the right-hand-side (R.H.S.) of 
([^ can no more be exact. 

In order to obtain an implementable scheme, a precise ODE integration for the 
characteristics and a quadrature rule can be used. We follow an approach very 
similar to [31] for variable coefficients. It consists in using Gaussian quadrature 
formula to approximate Q in regions where the involved functions are smooth. 

Remark 2.2. Indeed, in |31| . an other SLDG scheme is presented, but our form 
is equivalent to one form of SLDG as explained in [HI Proposition 4.5]. This may 
lead to different programming algorithms however. 
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2.1. Preliminaries. Let {xa}a=o,...k be the set of Gauss points in the interval 
(—1,1), with its corresponding weights {wa}Q=o,...fc {wa > 0), such that: 


^ k 

Vp € P 2 fc + 1 , / p{x)dx = ^ Wo,p{Xa). 

1 a=0 


In particular, we get on the interval li, 


p k 

yp&P 2 k+i, / p{x)dx = '^wIp{xI^) 

^ Q! = 0 


( 10 ) 


( 11 ) 


where x^^ := Xi + Xa^x = x^_i + |(1 + Xa)Ax and := ^Wa- 

To each set of Gauss points {x'^a}a=o,..,k in h, we can associate the corresponding 
Lagrange polynomials (dual basis) {^'^a}oc=o,..,k defined by 




l/.(a;) 


n 

0</3<fc 

/3^a 


X — X^ 
Xoi Xj^ 


For any u'^ GVk^ there exist coefficients ^ ^ ^ such that: 


M-l k 

u”(a;) = 

a—0 

In particular, the left-hand side of (|^ for (/j = becomes 

J u'^+\x)(pl,{x)dx = J u'^+^{x)ipl,{x)dx = uiywl,. 

n li 


( 12 ) 


(13) 


2.2. Definition of the scheme in the general case. Due to the discontinuities 
of u”, we separate the right-hand side of ([^ into several integral parts involving only 
regular functions: the R.H.S. of ([^ is approximated by the Gaussian quadrature 
rule on each sub-interval where u^{yx{—At)) is a regular function. 

For a given mesh cell li, we first consider the points {xi^q)i<q<p- (in finite num¬ 
ber) of the interval (xj_ i ), such that for 1 < g < pi, ?/a,; ^(—At) = x^.^_i 
for some £i^q G 1, and := a^i.Pi+i := ^i+i (see Figure j^. Then we apply 
the Gaussian quadrature rule on each interval Ji^q = {xi^q,Xi^q+i) and obtain the 
following quadrature rule, for any polynomial ip G Vk'- 


Xi,q + 1 


f u'^{yx{-^t))ip{x)dx = Y [ u'^{yx{-^t))p{x)dx (14) 

Pi k 

- (15) 

q—Q a^O 

with := ^(cci,,-!-! - Xi^q) and i* „ := Xi^q + 1(1 -|- Xa){xi^q+i - x^,,) = 

+ l ^ ^ ^i,q + l ~Xi,q ^ 


2 


2 
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^i.Q ^i-l/2 Xi^i Xi,2 + 



Definition of the scheme (operator Tb.At)' is the unique element of 14 satisfy¬ 


ing for all ip GVk, 

f u^+\x)ip{x)dx = (16) 

„■ — n — n __n 


M—1 Pi k 


2=0 q—0 a=0 


The scheme is made explicit by using formula (16) on each ip = if^^. The scheme 
equivalently defines an operator 7b,At such that 

n+1 ^ n 

u ^ = Tb.At U . 


In particular, if b is constant, then 7b,At = 71,At, and this is no more true if b is 
non-constant. 


Definition 2.1. For further analysis, let us introduce the following scalar product 
on 14 (where the index ”G” stands for the use of the Gaussian quadrature rule): 

M — l Pi k 

((P,4)g := ^ ^{S;q,a)- (17) 

2=0 q—0 a—0 

Then the scheme ( |16[ ) is equivalently defined by 

(u"+\v3) = {u^{y.{-At)),ip)G, V(^G Vfc. 


2.3. Stability and error estimate for constant drift coefficient. The weak 
form ([^ gives the stability of the scheme in the norm, at least in the case when 
b = const. Indeed, taking ip = in ([^ we get 

||^-+i||2, = („-(. _ 5At), r"+ 1) < ||zi-(. - &At)|4. 

where || • 1^2 denotes the norm on D and (., .) is the associated scalar product. 
Then, by the periodic boundary condition, ||u”(-— &At )|42 = ||m "|42 and therefore 

h"+'||L^<lk”||L- (18) 

This proof works only for b constant, however. 

For any w G L'^, we denote its projection on 14 by IIw, corresponding to the 
unique element of 14 such that 

||ui-nw ;|42 = inf \\w- fW^^. 


(19) 
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Remark 2.3. The function defined by Q corresponds to the projection of 

the function x —)■ u^{yx{—j^t)) on the space 14 ; 

u"+i=n(M"(y.(-At)), 

and, in the same way, we have = XIiio- 

We now recall a simple estimate for the projection on 14. 

Lemma 2.1 (Projection error). Let k > 0 and £ < k. If w G then 

||ir —ni (;|42 < 

where Ci{w) := 2 <!+i(\+i)! 

Proof. Let us write w = P + R where P is the element of 14 corresponding, on 
each interval 4, to the Taylor expansion of w centered at Xi and of degree £. We 
have \\w — nw||j ;,2 < ||ri; — P||l 2 = ||i?||L 2 < |np/^||i?|4oo. By the definition of R 
and usual Taylor estimates, we have ||i?||L°° < CiAx^~^^. □ 


Let v^{x) := v{tn,x) where v denotes the exact solution of Using the Te¬ 
stability of the projection, it is straightforward to show that — nu ”“''^|42 = 

||n(u”(- — bAt) — u"(- — 6 At )|42 < ||u” — u”| 42 , therefore we have 


By using Lemma 2.1 


this leads to the following known convergence result [5T] . 


Theorem 2.1. Let k > 0 and b be a constant. Assume the initial condition vg is 
1-periodic and in . Then, the following estimate holds: 

\\u'^-v^L^<\\u°-v°h^+CT^^, yn<N, (20) 

where the constant C depends only o/ |n| and k. 


Remark 2.4. By taking At = Ax this leads to an error estimate in 0{Ax^). 
However the examples (such as in Example 1) will show a numerical behavior in 
0{Ax^'^^) (as already remarked also in |31j 1. We refer to the recent work in |36j 
for more insight about this gap. 


2.4. Non-constant b: preliminary resnlts. For u G 14, the following approxi¬ 
mation result is central. It controls the error between the desired formula (|^ and 
the implementable scheme (16). 


Proposition 2.1 (Ganss quadrature errors). Let k > 0 and let b be of class 
Q 2 k +2 1-periodic. Then: 

(i) For all u G Vk, 


{u{y.{-At)),ip)G - {u{y.{-At)),ip) 


< C'AtAx^||M|42||i^|42, 


e Vk. 


where C >Q is a constant. In particular, we have, in the L^-norm: 

fb^AtU^ = u^+^ =Tb,AtU^+ 0{AtAx^\\u^L^), Vn>0. (21) 

(ii) For all u G Vk, for any fj in 1-periodic, 


{u{y.{-At)) - 'tp{y.{-At)),g))G - {u{y.{-Af)) - ip{y.{-At)),ip) 

<CAtAx2||u-4|42||v7|42+CMfc+i(4)Ax'=+i||^|42, Vy^GVk, (22) 
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where C >Q is a constant which depends only of k, and 

Mp{^) := max \\. 

0<r<p 

(Hi) For any regular ip G for any ip G Vk, 

{ip,(p)G = + 0(Mfc+i(-!/')Aa;'"+^||y)||i2). 

(iv) Furthermore, 3C > 0, for any ip G 1-periodic, 

WfbAti’-Tb.Atf^WL- < CMk+i{iP)Ax’^+\ 


( 23 ) 


(24) 


(25) 


Remark 2.5. Some assumptions can he weakened, for instance (i) and (ii) are 
still valid using that is in L°°, then in the error bounds (21) and (22) the 

AtAx^ term should be replaced by AtAx. However these bounds will be used in 
Section 3 and the form ( |21| ) and \22\ is preferred. Also, it is possible to prove 
that the error term 0{Mi^^i{ip)Ax^^) in (ii), {Hi) and (iv) can be improved to 
0(M2k+ii.i^)Ax^'^‘^) provided that ip G 


Proof of Proposition \2. 1\ 

Notice that the estimates of (i) and (Hi) are a consequence of (ii) (either by 
choosing ip = 0 to obtain (i), or by choosing At = 0 and m = 0 to obtain (Hi)). 
Then {iv) is deduced from {Hi) when applied to the regular function ipi{x) := 
ip{yx{-At)). 

The plan is first to prove {i), and then to generalize to {ii). Precise estimates for 
the (2fc + 2)nd derivative of cc —)■ u{yx{—At)) will be needed in order to estimate the 
error when using a Gaussian quadrature formula. In the following, we first bound 
the derivatives of a: —>■ yx{—t). 

Lemma 2.2. Assume that b G , for some k > 1, and 1-periodic. Let L := ||&'|jioo 
and let t gM.. Then x ^ y = yx{~t) *5 of class , 1-periodic, and 

( lly||L“(o,i) A 1 + I!^||l“ 1^1, 

I I||^?/||l~(o.i) < (26) 

\^and,ifk>2, || —j/||ioo(o,i) < C , Vg € {2,..., fc}, 

for some constant C > 0. In particular, all the previous derivatives are bounded on 
a fixed time interval t G [0, T]. 

Proof. We consider y as a function of the time t and of x. We can assume that 
X G [0,1] since we have yk+x{t) = k + yx{t) for all fc S Z and t,x gM.. We denote 
by = g|fcy the fc-th derivative of y with respect to x. 

Firstly, y{t,x) = x + b{y{s, x))ds and therefore, for x G (0,1), |y(t,a;)| < 
l + \\bh^\t\. 

For fc = 1 and 6 G we have ^g|y = b'{y)-^y and g|y(0) = 1, therefore 
||^y(t)| =exp(/p6'(y(s))ds) < 

For k > 2, we have 

|y('=) = (6'(y)y(i))('=-i) 
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Then we use a recursion argument for I = 1,... ,k. Let us assume that the spatial 
derivatives are bounded for 1 < ^ < fc — 1, with ||l=o(o,i) < . 

Then for fc > 2, the function / := i® bounded, with a 

bound of the form ||/(., t)||ioo(Q i) < for some constant C. By using 

the formula, for a given and fixed a;, 

yW(t) = e-^o 

Jo 

the fact that y^^\0) = 0 for fc > 2 and for s G [0, t] (or s € [t, 0] if t < 0): 

|e/>'(2/W)d«/(5)| < |5|fe-2gT|d 

we conclude that \y^'^\t)\ < C'|t|'=-i □ 


Lemma 2.3. Assume g > fc + 1, and u G Vk- On any interval J where u is regular, 
nq ^ 

P=1 

Proof. We first recall an expression for the g-th derivative of the composite function 
u{y), also known as ”Faa di Bruno’s formula” [T^ : 

P=1 ^ (aj), T,i 0‘j=p, J2i joij=q 


(y(l)/l!)ai . . . 


ai! • • • a,! 


Here the sum is limited to p < fc (instead of p < g) since u G Vk- 
Therefore, together with Lemma we obtain the bound 

k 


P=1 


(p)l 


L-^(y(J)) 


^iC.2 + -+0 

(“j). E]=i“j=p. E]=ii“3=9 


The case when 02 = • • • = Oq = 0 happens only if ai = p = q. Since g > fc + 1, 
and p < k, this case never occurs. Therefore, the power of At is at least 1, which 
concludes the proof. □ 


Proof of Proposition 2.1 i): Let e be the error term, defined by 

M—1 Pi k 

e-= / ;■ 

2—0 q—0 Q—0 

We have e = J2i YTq=o ^i,q where 


»1 M—l Pi k 

:= / u{yx(—At))p(x)dx -EEE<. u(ys,i^^^(-At))p(x\of). 

*^0 ,-_n ^_n n 


£i,q-= u(y,,(-At))p(x)dx-J2^l,au{yil,J-^^))pi^l.ci) (28) 

Q, = 0 

and with := {xi^q,Xi^q+i). 

Let u('y) be the function x -G u(yx(—At)). Since u(y) is regular on Jt^q 

for each fixed i, q G {0,...,pi}, and that the R.H.S. of (28) corresponds to the 
Gaussian quadrature rule on Ji^q, then we have in particular 

where := Xi^q+i - x^^q. 
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On the other hand, since (p G Vk, 

k 

r—0 

For all r G {0,..., k} we have 2k + 2 — r > k + 2 > k + 1, hence we can use 
Lemma 12.31 and obtain the bound 

k k 

r=0 p—1 

In particular, 


k k 


Pi 




r=0 p—1 i q—0 


By a scaling argument 0 na, and using that p G Vk for fixed k, we have, 
yO < r < k, 




Ax, 


c c 


(29) 


for some constant C, assuming also Axi^q < 1 (the idea is to use the fact that for 
polynomials of degree k, by using norm equivalences, H^cxj^q < C'||</5||i2('g 
for some constant C independent of p, and then to use a scaling argument from 
(0,1) to Ji^q to obtain the desired inequality). 

Denoting by | J| the length of any interval J, we have also 

< \yiJ.,q)\ < \J^,q\e^^\ L := \\b'h^, 


where \Ji,q\ = Axi^q. Hence, for r < k and p < k, 

A ™2fc+3 Il,„(r)|| ll„.(p)ll / /K ^2k+3\:^ 


i,q 


i,q ‘—^'^i.q 


< CAxlqY, 


Finally, by the Cauchy-Schwarz inequality. 


i^q ' i,q ^ ' i,q 


< 


l,q 

l^\\u\\l^. 


1/2 


since (J- ^ Ji^q is a covering of [0,1]. Hence we obtain 

^ \e^^q\ < CAtAa;^||(/)||L2||u||L2, 


which concludes the proof of (i). 

Proof of Proposition \2. l\ ii) : Let us write if = P + R where P G 14 is defined as 
the Taylor expansion of if on each Ji^q = {xi^q,Xi^q+i), around Xi^q. We consider 
the decomposition 


u{y.{-At)) - if{y.{-At)) = {u- P){y.{-At)) - R{y.{-At)) 


(30) 
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Then by Proposition i), for any ip G Vk, 

\{{u- P){y.{-At)),(p)o - {{u- P){y.{-At)),(p)\ < CAtAx'^Wu- P\\l 2 \\ip\\i^ 2 . 
Using the fact that ||-R||l 2 < C'||i?||Lc=o < we obtain the bound 

\{{u- P){y.{-At)),(p)G - {{u- P){y.{-Atj),ip)\ 

< + C'Mfc+i(i/))AtAa:''+^||(^||i2. (31) 

There remains to bound the error 


{R{y.{-At)),ip)G - (i?(y.(-At)),v?). 


This is easily bounded by C||i?||oo||<^|| l 2 = 0{Ax'^~^^\\<p\\l2). Combined with (30) 
and (31), we obtain the desired bound. □ 


2.5. Non-constant b: stability and error analysis. We now turn on the stabil¬ 
ity and convergence analysis. The following result shows the unconditional stability 
of the scheme, for any A: > 1. 


Proposition 2.2 (Stability). Let k > 0 and let b be Lipschitz continuous and 1- 
periodic. Then: 

(i) for any u G and u(x) := u{yx{—t)), it holds: 

11111^2 < ||u||i 2 , where L ■=\\b'\\L'^. (32) 

{ii) If furthermore b is of class there exists a constant Ci > 0 such that, 

Vm e Vk, 

\\fb,Atu\\L^ < Vu G Vfe. 

(in) In particular for the scheme = Tb,AtU^, 

||m”||l 2 < e'^i*"||u°||L 2 , Vn>0, 

where t„ = nAt. 


Proof, (i) We make use of the change of variable x —>■ z := yx{—t), with periodic 
boundary conditions for the integrands. Therefore we have x = yz(t) and 

r-t 


|^(t) = exp b'(yz(s))ds^ < 


We then obtain 


[ \u{yx{-t))\‘^dx = [ |it(z)P ^(t) dz < [ |u(z)pdz. 

Jn Jn (dz Jq 

{ii) By using (21), we have 

l|■^,AtM||L2 < \\u{y.{-At))\\L2 p CAtAx^\\u\\L2. 


Together with (32) we get a stability constant 

ef -h CAtAx^ < e^^\l + CAtAx^) < gf 
hence the desired result for any Ci > 0 such that Ci > \L CAx^. 


(33) 


□ 


We now state a first convergence result. It generalizes the error estimate of 
Theorem |2.1| established in the case when b is constant, to the non-constant case. 
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Theorem 2.2 (Convergence). Let k > 0. Assume the initial condition vq is 1- 
periodic and of class Let b he 1-periodic and of class There exist 

constants Ci > 0, C > 0 such that 

\\u--v^L^<e^^^(^\\v°-u°h2+CT^^y yn<N. (34) 


Proof of Theorem \2.^ By using the regularity of i;" 
have 


and Proposition 


2.1 


j) we 


= Tb,Atv^ = Tb,Atv^ + 0(Ax'=+'). (35) 

Because of the projection error = 0{Ax^~^^), then we obtain the 

following consistency estimate: 

= fb,Atv" + 0(Ax'=+i). (36) 


Therefore 


^n+l_^n+l ^ TbM^"- 

By the stability bound of Proposition |2.2 [ m), 

We conclude by induction. 


v^)+0{Ax^+^). 

'^\\l2+CAx'^+\ 


(37) 


□ 


2.6. Stability to perturbations. We conclude by a stability result with respect 
to the error of the position of the characteristics. 


Proposition 2.3. Let wi{x) := yx{—Af) and w^ix) := yx{—At) be some approx¬ 
imation of yx{—At) such that max|'u;i(ai) — x\ < CpAt for some constant Cq > 0. 

2=1,2 

Assume that ^ < K for some constant K > 0. Then for all u,ipG 14, it holds 


u{w 2 {x))(p{x)dx — / u{wi{x))(p{x)dx 


Ax 


11 “ 1111 11 J 




for some constant C > 0 independent of At, Ax. 


Proof. We first notice that \yx{—At)—x\ < c^At < Cq^Ax < qAx for some integer 
<7 > 1, as well as \yx{—At)—x\ < qAx. For a given interval /, let Iq := I+[—q, gJAx. 
It holds: 


|| u (' u ; 2 ) - u{wi)\\l2(i) 


< 

< 


|l~(/,)||w2 - wiIIl^ 




IWh^ilAu II A 1/2/ II II 11^2 

Cl ^^3/2 11^2 - williocAa: / < Ci||u|U2(j_^)- 


Ax 


for some constant ci > 0 (we have used a scaling argument as before). We remark 
that ||w||| 2 (/^) = J2j=-q,...,q ll'“lli 2 ( 7 + 5 Ax) where J = I+qAx is also another interval 
of same length as I. Hence I|ii|li 2(7 ) = (2g + l)||u||| 2 , and 


l{w2) - u{wi)\\l 2 < Ciyj2q + l||u||7,2 


\\W2 - Wi||l° 
Ax 


The result (381 follows by using a Cauchy-Schwarz inequality. 


□ 
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Corollary 2.4. We consider that an error is made in the computation of the char¬ 
acteristic yx(—^t), such that 


\y^{-At) - y^{-At)\ < e (39) 

for some constant C > 0 and e > 0. Then the error estimate of order 
Theorem 2.2 must be replaced by 


At 


in 


CT- 


At 


CT 


Ax At 


Sketch of proof. At each time step an error of order e = ||ui 2 ~ is made in 

the computation of the characteristics. By the previous Lemma this results in a 
supplementary error term of order Hence after At = -^ time steps the error 
coming from the computations of the integrals will be bounded by □ 


We remark that in practice, this approximation error is not seen in the numerical 
tests because the characteristics are computed using an analytical formula or a 
machine precision fixed point method when needed. A high-order approximation 
method would also lead to £ := CAt'^'^^ in (39) which can be made arbitrarily 
small in particular because we deal only with one-dimensional approximations of 
characteristics in the proposed method. 


3. Second-order PDEs 

This section deals with SLDG schemes for second-order PDEs. We will first 
deal with a simple diffusion problem with constant coefficients, for which specific 
schemes can be obtained, and then we consider the more general case of advection 
- diffusion problems with variable coefficients. 


3.1. Case of a diffusion equation with constant coefficient. We first consider 
a diffusion equation with a constant coefficient cr € K: 

„2 


vt - -^Vxx = 0, cc e D, t e (0, T), 
v(0,x) = vo(x), X G ft, 


(40) 

(41) 


and aim to construct simple schemes in this particular setting. Following Kushner 
and Dupuis |52], a first scheme, in semi-discrete form, is 


n+1 


(a;) = i ( u 


’’(x — a'/At) u"(a: -I- aV~At)^ = {x). 


(42) 


It is easy to see that, taking v"(x) := v(tn,x) where v is the solution of (40) and is 
assumed sufficiently regular, the following consistency error estimate holds: 

„,n-|-l cO ,,n 

11^^-^^IU^ = 0(At). 

The basic SLDG scheme (also called hereafter SLDG-1) is based on the weak for¬ 
mulation of (42). 

SLDG-1 scheme: Define recursively in 14 such that 

J u’^~^^{x)ip{x)dx = y i ^M"(a; — o-\/At)-f it"(a:-I-cr-s/At)^ (/?(a:) dx, V(/7 € 14 . 
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(The initialization of is done as before). The scheme will be also written in 
abstract form as follows: 

=5At(R"), 

where 

SAt ■= ■ 

Before doing the numerical analysis, our aim is first to improve the accuracy with 
respect to the time discretization. The technique proposed here is to use a convex 
combination of it, Sai, SAtSAt, ... It will work only for the constant coefficient case 
(cr constant). 

Using Taylor expansions, for u sufficiently regular, we have, for At small. 




At‘ 


SltSltU 


2 4 

U + + O(At^), 

4 

u + cr'^UxxAt + %-u^x'^ At^ + O(At^), 
o 


where u'x'^ denotes the g-th derivative of u w.r.t. x. 

On the other hand, if n" = v{tn,x) where v is the exact solution of Vt 
we have 


— u” + vtAt + —vttAt + O(At^) 


(43) 

(44) 


(45) 

(46) 


Now, looking for coefficients a, 6, c snch that av'^ + equal 

to up to 0{At^), using (43) and (44), we obtain the system 



c = 1 

1 


c = 

c _ 

3 “ 


(47) 


and we find that a = b = c = I. Therefore, a second-order scheme (for constant 
coefhcient) is now given by 
SLDG-2 scheme: 




= SltU^ ■■= + ^At^At^^”). 


(48) 


Remark 3.1. A variant of this scheme can be 


n+1 


1 




(49) 


This is in general slightly different from (48) because SAtSAt = n5'Atn5'At fnay 
differ from Nevertheless, the difference between the two will be of the 

order of the projection error 0(Aai^''’^) when applied to a regular data. 
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In order to obtain a third-order scheme, we can proceed in a similar way. First, 
we obtain the following expansions: 


5' 


0"^ . cr'* . 9 cr® 


At 


U = u- 


+ 0(A<4), 


= u-\- 

O 40 

‘S'it'S'Af'S'AtW = M+^cr^WxxAt+^t7‘^4'‘)At2-^^cr®4®^At3-tO(At^). 

Looking for coefficients a, b, c, d such that au" -I- 
is equal to i;”+^ up to O(At^), we find the system 


b 

b 

2 

24 

6 ! 


C + 

c -I- 

- + 

Q ' 


d = 
ld = 
id = 


Ac + ^d = 

45 ^ ^ 240 


1 

2 

1 

8 

48 


(50) 


and its solution 


(a,6,c,d) := —(13,21,9,2). 
45 


Thus, the following scheme is of 3rd-order in time: 

SLDG-3 scheme: 

= SltU^ := + ^^At^Atu" + ^SAtSAtSAtu'-^- 

45 15 5 45 

As in Remark 13.11 a variant of the scheme can be 


u 


n+l 


_ JJ _|_ ^ cO „,n I 1 cO cO „.n 2 




V45 


(51) 


Since we are using a convex combination of stable schemes (S'At, SAtSAt or 
SAtSAtSAt), the schemes SLDG-1, SLDG-2 and SLDG-3 are all stable in the 
norm. 

Remark 3.2. Up to 5th-order schemes - in time - can also be obtained (see m), 
using convex combinations of the form ^ii^AtY'^"^- 

We now state a convergence result for ([40|. 


Theorem 3.1. Let k > 0 and let a be a constant, and assume that the exact 
solution V of (40) has bounded derivative for q = max(/c -|- 2,2p -|- 2). IFe 
consider the SLDG-p schemes with p = 1,2 or 3. Then 




CT{ 


Aa;''+^ 

At 


+ AtP), yn<N. 


(52) 


Furthermore the same results hold for the variants (49),(51) forp = 2,3. 


Proof. We will consider the proof in the case of the SLDG-2 scheme, with p = 2, 
the other cases being similar. By using the regularity of the exact solution (^ 
and bounded), we have the following consistency estimate: 

z,"+i = aov^ + + 0(At3), 


(53) 
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where oq = oi = 02 = |, and the bound O(At^) is in the norm ||.||l 2 . Since 
= nS'^jn^+0(Aa:^'''^) for regular data ip, we have also 5'^jU" = n(S'^j)^u"+ 
0{Ax^~^^), and thus 

z;"+i = aou” + aiSAtv^ + a2^Ai5Atu" + 0{At^) + 0(Ax^+^). (54) 


By the definition of the scheme we have 




= ^a,(SAt)h 


(55) 


i=0 


We deduce, using the consistency estimate (53), 


|u»+i _^"+i 


|l= < llJ^a,(SAt)^(u^-v^)h-+CAt^ + CAx^+^ 

i<2 

< J^a,ll(SAt)'(u" - v")h2 + CAt^ + CAx^+^ 

i<2 

< Hu’^- v’^Hl2 +CAt^+ CAx'^+\ 


(since > 0 and = 1). The result follows by induction. 


□ 


3.2. Advection-diffusion with variable coefficients. We recall that for the 
following PDE: 


a(t, x)^ 


b(t, x)vx + r(t, x)v 


f{t,x), X €fi, t€ (0,T), 

(56) 


with n = IR and terminal condition v(T,x) := w{T,x), introducing a probability 
space (Q,F,P) with a filtration {Ftjoo, and a one-dimensional Brownian motion 
(Wt)t>o, and the solution Xs = of the stochastic differential equation 


dXs = h(s,Xs)ds-\-<7(s,Xs)dWs, s>t, 

Xt = X, 


and if u is a regular solution of the PDE (56) on (t,T) (assuming that the partial 
derivatives dtv and exist and are continuous) then the following equivalent 
expectation, or ” Feynman-Kac” formula, holds: 


j{t, x) = E 


^-J'P'r{e,Xa)d0 


fT 

(T, AV) + e- xl’^) ds I At 


(57) 


To simplify, we shall focus here on the case when b and a do not depend of time, 
and r is constant. We consider the forward PDE: 

Ut - ^Uxx - b{x)ux + ru = f{t,x), xen,te{0,T). (58) 

In that case the Feynman-Kac formula gives, with h = At, T = t + h and u^(x) := 

u{tn,x): 


u'^+^ix) = E 






At 


■ ui(/l, x) 


(59) 
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with 


wih^ x) := E 


e-^^f{tr, + h-s,X°nds I Tt 


(60) 


/ -.2 

Let Aw := + b{x)wx — rw. The term w{h,x) is also the solution at time 

s = h oi the linear problem wt{s,x) = {Aw){s,x) + f{s,x) with initial condition 
w{0,x) = 0, and with f{s,x) := /(t„ + s,x). Assuming that the source term 
/ is regular and that we can use its derivatives, we can approximate it with an 
error by using a Taylor expansion: w{h,x) ~ Z)i=i (where 

Wjt denotes the j-th derivative with respect to time). In particular, wt{0,x) = 
f{0,x) = f{tn,x), and wu = {Aw + /)* = Awt + /t = A{Aw + /) + ft, so 
Wtti0,x) = {Af){tn,x) + ft{tn,x). Hencc in order to devise a second-order scheme 
we approximate (601 by 


){h,x) = hf{tn,x) + Y{-^f + ft){tn,x) + 0{h^). 


(61) 


The modihcation of the scheme is obtained, therefore, by adding at each time step 
the following correction term atGGauss quadrature points 


hf{tr, 




< x) + ~^{^f + /t)(^ii 


r). 


(62) 


For the approximation of the expectation in (59), we aim to use a higher-order 
semi-discrete approximation also called ’’weak Taylor approximations” in the sto¬ 
chastic setting, see in particular Kloeden and Platen Chapter 15]. General 
semi-discrete (and fully-discrete) approximations can be found in [ 22 ) . 

We will focus on Hrst- and second-order weak Taylor approximations. Some of 
these approximation may use the derivatives of b and a (Milstein [25], Talay |37|, 
Pardoux and Talay [28]). In our case we shall use a derivative-free formula of 
Platen |5D| (explicit second- and third-order derivative-free formula can be found 
in Kloeden and Platen [20], as well as multidimensional extensions). 

Let us denote b — b{x), a = (t(x) as well as 7 ^^ = ^\^{x)-. 

7At(*) X-\-b{x)Xt-\-q(7{x)'J^i. (63) 

Our SLDG-1 scheme, corresponding to a first-order (weak Euler scheme), is 
defined by 




= S 


(i)„ 

At^ 


:= n 


\y\t{-)) 


(64) 


'g=±l 

with weights a_i = ai = ^ and characteristics = it,- 

Our SLDG-2 scheme, corresponding to the second-order Platen’s scheme, is de¬ 
fined by 


1 - 1-1 — 


At ' 


:= n 


Y «9«"(yAt(-)) 


(65) 


with weights a-i = ai = | and ao = | and characteristics = y\{x) defined by: 


yl{x) = a:-F^( 6 ( 7 , 


f') 


b)h 


(o-(7i^) +o-(7/i^) + 2cr)v^ q+{<7{-il) - a{'^^^)){M^ - 1) 


( 66 ) 

Vh. 
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Remark 3.3. In the constant coefficient case cr(a:) = a, the scheme becomes 

y^n+i ^ := — aVsKt) + + crVSAt) j (67) 

Remark 3.4. Higher-order weak Taylor schemes can he found in [20] and could be 
used with DG to devise fully discrete schemes in the same way. 

The above SLDG-1/2 schemes are no more exactly implementable because b{x) 
and cr(x) are not constant. So, as in the advection case, we consider the use of a 
Gaussian quadrature rule on each interval of regularity of the data. 

Remark 3.5. Notice that if h is small enough such that 

II+ veer'll< 1, (68) 

then for each q = ±1 the function x —)■ 'j'^(x) is a one-to-one and onto function. 
Furthermore, its inverse can be easily and rapidly computed by using a fixed point 
method or Newton’s algorithm. Details are left to the reader. 

In the same way, for h small enough such that, for instance, 

^I1^'I|l“ + 3v^||ct'||l=o < 1, (69) 

then X — )■ y®(x) as defined in ( [6^ is one-to-one and onto function. 

SLDG-1 scheme (fully discrete): For each given p = ±1, we consider a partition 
of li into intervals such that all are subintervals of some Ij. We then 

define Gauss points and the bilinear product (a, b)cv in a similar way as in ( [T7| ), 

that is, using the Gaussian quadrature rule on each Hence we define 
in 14 such that 


e 14. 


(70) 


ri=± 


Formula (701 involves two different quadrature rules, because the discontinuity 
points of M"(y“'"(x)) and m”(j/“(x)) are not the same. It differs from the definition 
of which satisfies 


V?) = ^ e Vfc. 


(71) 


T] = ± 


SLDG-2 scheme (fully discrete): In a similar way, we define 5'Xt m” in Vk by: 


(>S'a1m”,(/?) = X «'7(w"(yAt)>‘7’)G’7, VtpeVfc. 


(72) 


3.3. Stability and convergence. We first state some useful estimates for t he 
operators SAt € S'^t}. The proof is similar to the one of Proposition 


2.1 


Proposition 3.1. Let k > 0 and let a be of class and 1-periodic. Then: 

(i) there exists a constant C > 0 such that, for any y‘^^, for all m S 14) 


'MyAt)^T)Gi - {uiyAt)^T) 


<CVMAx^\\uh4T\\L^ 
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In particular, for any u G Vk, 

SAtu = SAtu + 0[VAtAx'^\\u\\L2). (73) 

(ii) For all u G Vk, for any fj in 1-periodic, 

SAtiu -Ip) = SAtiu -Ip) + 0 {VMAx'^\\u - V’IIlO + 0{Mk+i{ip)Ax^~^^), (74) 

where C >Q is a constant. 

(Hi) For any regular ip G , 1-periodic, we have in the Lf norm 

SAti’ = SAti^ + 0{Mk+i{'iP)Ax^+^). (75) 


We now establish stability properties. 


Proposition 3.2. Let k>0, and assume that h is small enough in order that (68) 
(resp. holds. 

(i) (Stability with exact integration as in (64).^ For any m S 14, 

WSaMIl^ < (1 + C'At)|jM||i2, 
where C > Q is a constant. 


(ii) (Stability with Gaussian quadrature rule as in (70).} For any uGVk, 

WSaMIl^ < (1 + CAt + cVMAx^)\\u\\l2. 

(in) In particular the fully discrete schemes SLDG-1 and -2 are stable under 
the ’’weak” GFL condition 


Ax'^ < XAt, for some A > 0. 


(76) 


Proof, (i) By making use of the convexity ol x ^ x^, the change of variable formula 
X —>■ y\^(x) (and denoting also z -G x\^(z) the inverse function of we have 

2 


ISaM 


L2 


aqu(x + b(x)At -|- qa(x)\/At) 


< 


E 

Q 

E 


u(x + b(x)At + qa{x)'/At) 


dx 


dx 


1 + b’(x'i(z))At qa'(x‘i(z))^fAf 


\u(z)(^dz. 


Then we remark that x\^(z) = x + 0(-\/At), so 1 + b'{x^{z))At-\-qa'(x‘^{z))'/Ai = 
l-\-qa'(x)'/At-\-0(At), and for At small enough 0 < (l-\-b'(x‘‘(z))At-\-q(T'(x'^(z))\fAty 
1 — qa'{x)'/Ai + CAt for some constant C > 0. Hence 


||5'Atw|| 


L2 S 


[ 0 ^ 9(1 ~ q(T'(x)'/At + CAt)\u(z)(^dz 

< (1 + CAt) J \u(z)ydz 


where we have used that Y^aq = 1, and Yq = 0- The desired result follows. 
(ii) This is a consequence of (i) and of the bound (73) of Proposition 3.1 □ 


The convergence result for the approximation of (58) is the following. 
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Theorem 3.2. Let k > 0 and let a be a 1-periodic function, of class We 

consider the schemes SLDG-p for p = 1,2 (implementable version). 

Assume the exact solution v has a bounded derivative for q = max(2p+2, fc + 1), 

and that the weak CFL condition (76) is satisfied, then 


LiT 


<e 

for some constant Li > 0. 




■CT{ 


At 


+ AfP) , yn<N, (77) 


In particular for At = XAx for any A > 0, and k = p G {1,2}, the SLDG-p 
schemes are fully discrete schemes and of order 0{AxP). 

Proof of Theorem \3.2\ We first consider the SLDG-1 scheme = S/^t'uA. By 
making use of the consistency error estimate, we have 

^"+1 = -k 0{At^) + 0(Aa;'=+i) = 5'Atu" -f 0{At^) + 0{Ax'^+^). (78) 

Furthermore, by proposition |3.1} iM), 

- 5Atu"|U2 < CMu+i{v^)Ax'^+\ (79) 

Hence 

v'^+^ = SAtv^ + 0(Af2) + 0(Ax'=+i), (80) 

and by difference with the scheme = S'Atu": 

^ \\SAtu^-SAtv^L-+C{At^ + Ax'^+^) (81) 

< +C(At2 +Ax'^+i), (82) 

for some constant C > 0, where we have made use of the stability estimate for S'ag 
T herefore we obtain the desired error bound. 

For the SLDG-2 scheme, the estimates are similar, using the fact Platen’s scheme 
is second-order to get the consistency estimate u”-|-0(At^)-|-0(Aa;^'''^). 

The conclusion follows. □ 


4. Extension to two-dimensional PDEs and splitting strategies 

4.1. First-order PDEs - two-dimensional case. We aim to extend the previ¬ 
ous scheme to treat two-dimensional PDEs, by using splitting strategies and one¬ 
dimensional solvers of the previous section for advection in the direction of the 
coordinate axes. 

Let H be a square box domain D = [xi^rnin,Xi,max] x [x 2 , min, x 2 , max] with pe¬ 
riodic boundary conditions. Let us consider a spatial discretization of LI into cells 
lij := li X Jj where li (resp. Jj) is a cell discretization of [xi^min,Xi^max] (resp. 
[x 2 ,min,X 2 ,max]) as in the one-dimensional case using Mi (resp. M 2 ) points. We 
define the corresponding space of 2d discontinuous Galerkin elements by using the 
Qk basis (v G Qk if v{x) = SjjXfe %a:}a;^): 

:= G L\n,R), G Qk, V(z,j)|. (83) 

We consider the case of 

Ut + bi{xi,X2)ua;^+b2(xi,X2)ua;2=0, {xi,X2)gLI. (84) 
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The idea, already proposed in [34] or |9] is to split the equation into 

ut+ bi{xi,X2)uxi = 0, {xi,X2)€Q (85) 

and 

Ut +b2ixi,X2)Ua:^ = 0, {xi,X2)€^- (86) 

Let the corresponding characteristics defined by : 

• for g = 1: ,^^)(t) = (yi{t),X 2 ) where 

yi{t) is the solution of yi{t) = bi{yi{t),X 2 ) with ?/i(0) = xi, 

• for g = 2: = {xi,y 2 {t)) where 

y 2 {t) is the solution of y 2 (t) = b 2 {xi,y 2 {t)) with 1 / 2 ( 0 ) = X 2 - 

Let be the corresponding exact evolution operator in the direction of Xq. The 
exact solution of (85), with g = 1 (resp. (86), with g = 2) satisfies 

v^+'^{xi,X2) = i;”(X‘(^^^^^)(-At)) = S%^{v^){xi,X2). 

We define the discrete evolution operator for (85), denoted so that for 

each fixed Gauss points X 2 = x]^ the one-dimensional scheme is used for the evo¬ 
lution in the direction cci. We define in the same way the operator for the 

approximation of (86). 


Remark 4.1. In the case of (85) we do not try to compute precisely the 2d integrals 
[ u'^{X'l^^.^^^{-At))(pi{xi)ip2ix2)dxidx2, (87) 

J hx Jj 

where Lpi and (f 2 are polynomial basis functions. The discontinuities of the inte¬ 
grand are no longer well localized and it would not be possible to obtain easily an 
accurate approximation for (87). Rather, the discrete scheme computes a high-order 
approximation of the following integrals on a full band [0,1] x Jj 


J[o,i\xjj 

and this is all what is needed. 


u''{Xl,,,^,,,^^{-At))ipi{xi)<p2{x2)dxidx2, 


( 88 ) 


Now, the results of Section 2, in particular Propositions 2.1 and |2.2[ can be 
extended to the operators ^^,q = 1,2. The difference is now that the consistency 
estimates are typically as follows, for g = 1, 2: 

< CA^AxI+^Ml^, V<^ g vI:^\ 

and 

Let furthermore £t be the evolution operator for the initial advection prob¬ 
lem (84). In the case when b = (&i, ^ 2 ) is constant we have 

= ^At^At 

and we can therefore approximate the exact evolution SaiV^ by 7),^ AtT~bi At'^^ with 
no error coming from the splitting. 

In the following, when there is no ambiguity, we furthermore denote 


n~Q _ q~Q 

'At ~ 'ba^At 


9 = 1,2. 
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In the case when b — (&i, 62 ) is non-constant, we recall the following approxima¬ 
tions of the exponential for A and B matrices and for small At: 

g(A+s)At ^ gSAtgAAt (Trotter Spitting), (89) 

g(A+B)At _ gS^gAAtgS^ _l_ (Strang’s spitting). (90) 


leading us to consider the following splitting approximations 

TbAt - TltTit (Trotter) (91) 

TbAt-TLTljL (Strang) (92) 

2 2 

of expected consistency error 0(At) and 0{At'^) respectively!^ These last two 
splitting schemes are similar to the ones used in Ell- 

Following [33], we shall also consider a 3rd-order splitting scheme of Ruth |55] . 
a 4th-order splitting scheme of Forest m (see also Forest and Ruth HH), as well 
as a 6 th-order splitting of Yoshida jUj). 

Ruth’s 3rd-order splitting: 


%At 


'j~2 'j~l 'y'2 -y"! 'y'2 

' ciAt ’ d\At ‘ C2At ’ d2At ' c^At ’ d^At"! 


(93) 


with 


Cl = 7/24, C 2 = 3/4, C 3 = —1/24 and di = 2/3, ^2 = —2/3, = 1. 

Forest’s 4th-order splitting: 


with 


IbAt — '72At '(^^+^, 3 )^ '72At '(^^+^ 2 )^ '72At ' 71 ^’ 


1 , 21/3 


(94) 


Yoshida’s 6 th-order splitting: 

o- ^ ^Ath ‘nrAth ^Ath 
'bAt — ' yiAt 'y2At 'y^Ati 

where denotes the previous Forest’s 4th-order approximation method, 


(95) 


1 


2/1 := 


2 - 21/5 


and j /2 := — 


21/5 

2-21/5' 


Remark 4.2. Stability in the LS-norm is then easily obtained. Indeed, we have 
the L^-stability of the one-directional advection operators that is, for variable 
coefficients 

\\rlM\L^<e^^^\\uh2 (96) 

for some constant c. Then, for instance for the Trotter splitting, we/lare 
g2cAt||y||^2^ w/iic/i gives the if stability result 

(97) 

In the same way any finite product of operators of the form (or any convex 

combination of such products) would lead to stable schemes. 


iDenoting r = T/N for AI > 1, and (J > 0, if linear operators At and Bt on a normed 
vector space satisfy At = Bt -\- 0(r'i+^), with ||A”||, ||B"|| < C for all 0 < n < N, then 
Af = Bfi + 0 {t1). 
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Hence the results of Section 2 can be extended: for a = 1,2,3,4 and 6 corre¬ 
sponding to the splittings (91), (92), (93), (94) and (95) respectively, for regular 

(98) 


solutions, the one time step error will be of order 

0(At“+i) + 0(Ax'=+i), 

and the convergence error bound after N time steps will be of order 


0(At° 


^ At ’ 


(99) 


4.2. Second-order PDEs - two-dimensional case. We consider the case of 
Ut — ^Tr{a{x)a{x)^D^u) + b{x) ■ Vu = f{t, x), x £ fl, t £ (0, T) 

( 100 ) 

(with initial condition u{0,x) = Uo{x)), where cr(a;) € and Tr{A) denotes the 

trace of the matrix A. 

We introduce the following decomposition into the direction of diffusions repre¬ 
sented by the column vectors of the matrix a (similar decompositions have been 
used by Kushner and Dupuis [22] , Menaldi [24] , Camilli and Falcone [4] , Debrabant 
and Jakobsen m , etc.): 


aa 


= cTqa^, where aq := 


9=1 


0'2,q 


Setting Bi = 


and i?2 = I ^ I, we write (]100[) as follows: 


Ut 


+ (-iTr{aqa'^D'^u) + Bq-\/v\=f{t,x). 

q^l,2 ^ ^ 


( 101 ) 


Let us first consider the one-directional problem (one direction of diffusion): 


1 . 


Ut - -Tr{aqa^D u) + Bq 


Vm = 0. (102) 

For this subproblem we consider weak Taylor schemes exactly as for the one¬ 
dimensional SLDG-1 and SLDG-2 schemes (63)-(64) and (65)-(66|. Indeed these 
approximations are known to be also of order 1 and 2 in time for (102) in any 
dimension [20] . 

It remains to give the definition of a scheme, of sufficient order, for the approxi¬ 
mation in two dimensions for terms of the form 


n(ii"(2/L(-))) 

where H is the projection on and y\i{x) is now a vector of 1 


(103) 


Remark 4.3. In view of the definition of the characteristics (63) or (66), a typical 
problem is to compute accurately the projection on of a function of the form 


{xi,X2) -£ u'^ifl{h,Xi,X2),f2{h,Xi,X2)), (104) 

with h = '/At, where f\ and f 2 are regular functions with known expressions, and 
such that 


fi{0,Xi,X2) = Xi and /2(0,Xi,X2) = X2. 


( 105 ) 
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A high-order approximation of the term (103), or (104) in the general case can he 


obtained by using the PDE satisfied by v{s, Xi, X 2 ) := w”(/i(s, Xi,X 2 ), / 2 (s, Xi,X 2 ))- 
More precisely, assuming that it" is a regular function, we observe that dgV = 
{dsf, Vit"(/i,/ 2 )) and Wv = / 2 ) (where Df := (|^) and Vii = 

Therefore dgV = {dsf, {Df'^)~^Vv) = {Df~^dsf, Vv) and v is solution of 
the'PDE 

dsV - {Df~^dsf, Vv) =0, s > 0, (106a) 

i>(0, a;i, a;2) = u'^{xi,X 2 ) (106b) 

(the matrix inverse Df{s,xi,X 2 )~^ is well defined for small s > 0 since by the 
assumptions (105) we have Df{ 0 ,xi,X 2 ) = Id). Then we have a problem of the 
form (84) and we can apply the splitting approaches of Section 4 .I to obtain a 


high-order approximation of (104) on a DG basis. 


Remark 4.4. In the present work we will consider only numerical examples involv¬ 
ing terms of the formIlvT{fi{h,xi,X 2 ),X 2 ) orIiu^{xi, f 2 {h,xi,X 2 )) (i.e. f 2 {h,xi,X 2 
X 2 , or fi{h,xi,X 2 ) =xi), oroftheformIlu"'{fi{h,xi),f 2 {h,X 2 )) with regular func¬ 
tions fi and f 2 and h = VAt. For such cases, the one-dimensional discretization 
can be extended to two dimensions by straightforward splitting. 


Finally, for the general case of (100), we define the scheme by using Strang’s 


splitting of the one time-step evolution operators for (102) and by adding the cor¬ 


rection (62) for the source term. 


5. Numerical examples 


The first three examples are devoted to advection problems, while the other 
examples concern second-order equations. 

We recall that N is the number of time steps (and At = T/N), and M is the 
number of spatial mesh points in the one-dimensional case (resp. Mi, M 2 for two- 
dimensional cases). 

Unless otherwise specified, the characteristics are one-dimensional and are always 
computed exactly (see added sentence in Section 5 before the first example). 

Computations were performed on a DELL Latitude E6220, Intel Core i5, 2.50GHz, 
4GO RAM, with Linux OS, 32-bit, using GNU C4—h. 


Example 1. We consider an advection equation with non-constant advection term 
Vt + h{x)vx = Q, a; G (0,1), t G (0,T), (107) 

d( 0, a;) = sm(27ra;), a: G (0,1), (108) 

and 


b{x) := Co -\- Cl sin(27rx), with Co = 1 and Ci := 0.8 (109) 


together with periodic boundary conditions on (0,1). The exact solution is given 
by v{t,x) = sin(27ri/a;(—t)), where 


, , 1 / / ,tan(7ra;) -I- r \\ 

yx[—t) = —atani — r -|- tani atan(- ) — CoTrat I I 


with r := & and a := Cl — r^. 

Co 

The results are given in Table for At ~ Aa; with fixed CFL= 1.8 and terminal 
time T = 1.3. (Here the CFL corresponds to ||6||oo^-) The numerical error 













SLDG SCHEME FOR FIRST- AND SECOND-ORDER PDES 


25 


behaves approximatively one order better than the expected one when At = XAx, 
that is of the order of ) = 0(Aa:^). Super-convergence results can be 

explained in some cases for other DG methods 9) . 



error 

k = 

1 

k = 

2 

k = 

3 

k = 

4 

M 

N 

error 

order 

error 

order 

error 

order 

error 

order 

10 

10 

1.95E-01 

- 

3.45E-02 

- 

1.45E-02 

- 

7.83E-03 

- 

20 

20 

2.67E-02 

1.93 

6.06E-03 

2.50 

1.38E-03 

3.39 

2.33E-04 

5.07 

40 

40 

7.80E-03 

1.77 

6.39E-04 

3.24 

3.22E-05 

5.42 

4.31E-06 

5.75 

80 

80 

1.47E-03 

2.40 

3.62E-05 

4.13 

1.52E-06 

4.40 

7.74E-08 

5.80 

160 

160 

2.27E-04 

2.69 

3.31E-06 

3.45 

7.13E-08 

4.41 

2.48E-09 

4.96 

320 

320 

3.92E-05 

2.53 

4.03E-07 

3.04 

3.92E-09 

4.18 

8.03E-11 

4.95 


Table 1. (Example 1) non-constant advection, At ^ Aa; and 
CFL= 1.8, T= 1.3. 


Example 2 (2D advection with non-constant coefficients). We consider the 
following rotation example of a ” bump”: 

Ut + 27 r(—ai2,a^i) • Vu = 0 , x = {xi,X2) G ft, t G ( 0 ,T), 

u{0,x) = 1 - 

with ft := (—2,2)^, tq = 0.25 and terminal time T = 0.9. Since b{xi,X 2 ) = 
2 tt{—X2,xi) is non-constant, Trotter’s splitting is no longer exact. 

In Table we test and compare the splitting algorithms as described in sub¬ 
section 2.3, from order 2 to 6 (Strang’s splitting. Forest’s 4th-order splitting and 
Yoshida’s 6th-order splittings, tested with k = 2,4, and k = 6 respectively), using 
Ml = M 2 = M spatial mesh points. Trotter’s splitting error, not represented in 
Table ([^, is of order 1. We have avoided taking the particular case of T = 1 (full 
turn) because it gives better numerical results but prevents proper understanding 
of the order of the method. 

In this example, the initial datum is sufficiently close to 1 outside a ball of radius 
1.5, so that the error coming from the boundary treatment is negligible. 



error 

Strang (with k 

= 2) 

Forest (with k 

= 4) 

Yoshida (with k 

= 6) 

N 

M 

error 

order 

cpu(s) 

error 

order 

cpu(s) 

error 

order 

cpu(s) 

10 

10 

2.91E-01 

- 

0.004 

1.66E-01 

- 

0.01 

1.8IE-02 

- 

0.07 

20 

20 

6.62E-02 

2.13 

0.012 

i.oiE-02 

4.04 

0.03 

2.45E-04 

6.21 

0.26 

40 

40 

1.60E-02 

2.05 

0.032 

6.24E-04 

4.01 

0.22 

3.64E-06 

6.07 

1.65 

80 

80 

3.99E-03 

2.01 

0.272 

3.89E-05 

4.00 

2.04 

5.61E-08 

6.02 

15.06 

160 

160 

9.96E-04 

2.00 

2.844 

2.43E-06 

4.00 

18.25 

1.03E-09 

5.77 

120.98 


Table 2. (Example 2), 2D rotation, L? errors at time T = 0.9, 
using M X M grid points and splittings of order 2,4 and 6. 


Example 3 (2D deformation with non-constant coefficients) In this exam¬ 
ple, close to the one in for instance Qiu and Shu ISl Example 5], the advection 
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term is non-constant 


Ux + (^ 9 {t) cos(y) sin(x)^'Uj^ = 0, 
ix,y) e 17, t e (0,T), 

with 17 := (—2,2)^, T = 1 and same initial datum as in Example 4. Here we 
furthermore consider g{t) := 1 for 7 S [0, and then g{t) := —1 for t g]|^,T], so 
that the exact solution after time T is m(T, cc, y) = uo{x^ y). 

In Table we test and compare the splitting algorithms of orders 2,4 and 6 
(Strang’s, Forest’s and Yoshida’s splittings), using polynomials of degree fc = 2, 4 
and 6 respectively. The cpu times are also given in seconds. 




error 

Strang (with k 

= 2) 

Forest (with k 

= 4) 

Yoshida (with k 

= 6) 

N 

M 

error 

order 

cpu(s) 

error 

order 

cpu(s) 

error 

order 

cpu(s) 

10 

10 

1.28E-01 

- 

0.005 

7.82E-03 

- 

0.08 

7.70E-04 

- 

0.85 

20 

20 

1.45E-02 

3.14 

0.034 

2.78E-04 

4.81 

0.36 

6.60E-06 

6.87 

3.65 

40 

40 

1.44E-03 

3.33 

0.104 

9.06E-06 

4.94 

1.58 

3.32E-08 

7.64 

16.20 

80 

80 

1.66E-04 

3.12 

0.620 

3.30E-07 

4.78 

7.73 

2.71E-10 

6.94 

140.41 


Table 3. (Example 3) 2D deformation, errors at time T = 1, 
using M X M grid points and splittings of order 2,4 and 6. 


Example 4 (ID convection diffusion). Now, we consider the diffusion equation 

vt - ^cr'^Vxx + bvx = 0, Va: G 17, 7 e (0, T) (110) 

u(0, x) = cos(27ra;)-I-^ cos(47ra;), cc G 17 (HI) 

together with periodic boundary conditions on 17 = (0,1), with constants a = 0.1, 
& = 0.3, and T = 0.2. The exact solution is given by 

v{t,x) = ^ Cfc exp(—2(T^fc^7r^7) cos(2A:7r(x — &7)), 

fc=l,2 


with Cl = 1 and C2 = 5. 

Since the operators and bdx commute, we use the simple scheme 

= S^^tTbAtU^. 

In Table we study the orders of the SLDG-RKp schemes when A7 ~ Ax and 
p G {1,2,3}. The orders are as expected. 

We also give in Table the errors when taking larger time steps (A7 ^ Ax), 
still showing good behavior, while the ratio ^ varies from 0.40 to 6.40. 

We have numerically also tested the case when 6 = 0 (pure diffusion); the nu¬ 
merical results are very close to the present case. 
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error 

SLDG-RKl (Pi) 

SLDG-RK2 (P 2 ) 

SLDG-RK3 (P 3 ) 

M 

N 

error 

order 

error 

order 

error 

order 

10 

10 

9.94E-03 

- 

1.37E-03 

- 

8.66E-05 

- 

20 

20 

1.39E-03 

2.84 

1.08E-04 

3.67 

3.70E-06 

4.55 

40 

40 

2.93E-04 

2.25 

3.63E-06 

4.90 

1.03E-07 

5.17 

80 

80 

8.02E-05 

1.87 

6.28E-07 

2.53 

9.81E-09 

3.39 

160 

160 

2.35E-05 

1.77 

9.72E-08 

2.69 

7.00E-10 

3.81 

320 

320 

8.22E-06 

1.52 

2.60E-08 

1.90 

5.79E-11 

3.60 

640 

640 

4.06E-06 

1.02 

6.17E-09 

2.08 

5.81E-12 

3.32 

Table 4. 
At ~ Ax. 

Example 4 (ID diffusion), SLDG-RKp schemes with 


error 

SLDG-RKl (Pi) 

SLDG-RK2 (P 2 ) 

SLDG-RK3 (P 3 ) 

M 

N 

error 

error 

error 

20 

10 

1.37E-03 

4.34E-05 

1.79E-06 

40 

15 

5.13E-04 

6.87E-06 

1.41E-07 

80 

20 

1.39E-04 

1.40E-06 

l.llE-08 

160 

25 

1.05E-04 

1.83E-07 

5.20E-10 

320 

30 

8.49E-05 

6.14E-08 

3.09E-11 

640 

35 

7.26E-05 

4.35E-08 

1.15E-11 

1280 

40 

6.35E-05 

3.31E-08 

7.02E-12 


Table 5. Example 4 (ID diffusion), SLDG-RKp with large time 
steps At Ax. 


Example 5 (ID Black and Scholes and boundary conditions) This example 
deals with the one-dimensional Black-Scholes (B&S) PDE for the pricing of a Euro¬ 
pean put option with one asset [38] . After a change of variable in logarithmic coor- 
dinatesj^the equation for the European put option becomes on ft := (xmin,Xmax)'- 

{ ut —\(j'^Uxx + bux+ru = Q, X e ft, t e {0,T), 
u(0, a;) = uo(a^) = AT max(l — e“, 0) x G Q, /-ini 

u{t, x) = ui(t) = Ke~^* - Ke^ t € (0, T), a; < Xmin, 
u{t,X^ - Uxit') =0 t € (0,T), X ^ Xyyiax, 

with b := —(r — ^cr^) and where Xmin < 0 and Xmax > 0, and we have imposed 
boundary conditions outside of Q. Numerically, the initial datum exhibits singular 
behavior at a; = 0 (as it is only Lipschitz regular). 

For this PDE the scheme reads 

ra+l _ -rAtqa- q-b n 

^ The classical B&S PDE for the put option reads 

vt — -a'^s^Vss — bsvs + rv = 0, s £ (0, oo), t S (0, T), 

(where b = r— with initial condition n(0, s) = ¥>(s) = max(iC — s, 0). Then using the change 

of variable x = log(s/iC) and u(t,x) := v{t,s), we obtain the PDE | [IT^ on fc G M. 
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The following financial parameters are used: K = 100 (strike price), r = 0.10 
(interest rate), a = 0.2 (volatility), and T = 0.25 (maturity). Since the interesting 
part of the solution lies in a neighborhood of cc = 0 (notice that has a singularity 
at a; = 0), for the computational domain we consider 

^ 2 ). 

In principle the PDE should be considered with \xmin\, |a^max| ^ 1, but here it can 
be numerically observed that the solution doesn’t really change for \xmin\, \xmax\ > 2 
Results are reported in Table for the errors, where At is chosen of the 
same order as Ax, and the SLDG-RKl SLDG-RK2 and SLDG-RK3 schemes are 
compared, together with a P 4 polynomial basis {k = A). We used a P 4 basis so that 
the error from the spatial approximation is in principle negligible with respect to 
the time discretisation error. We numerically observe the expected order 1 (resp. 
2) for the SLDG-RKl (resp. SLDG-RK2) scheme, and approximatly order 3 for 
the SLDG-RK3 scheme (of expected theoretical order 3). 

Remark 5.1 (Boundary treatment). For semi-Lagrangian schemes, the knowledge 
of u{t, x) for X < Xmin or X > Xmax Can be used if it is available. Here, ’’out- 
of-bound” values are needed for computing and for i;" = 

In particular, the values u^{x -\- ka^fAt — bAt) for |fc| < 3 are used when 
y := X -\- kaVAt — bAt lies outside of (xmimXmax)- In that case, we simply directly 
use the ’’out-of-bounds” values Uf,(tn,y) wheny < Xmin or Ur(tn,y) wheny > Xmax- 
It is clear that this will not work for a general PDE posed on a given domain with 
given boundary conditions. (See however [T] for an example of a semi-Lagrangian 
.scheme applied to a PDE with Neuman boundary conditions.) 


error 

SLDG-RKl 

SLDG-RK2 

SLDG-RK3 

M 

N 

error 

order 

cpu(s) 

error 

order 

cpu(s) 

error 

order 

cpu(s) 

10 

10 

6.30E-02 

- 

0.001 

3.84E-02 

- 

0.001 

4.17E-02 

- 

0.004 

20 

20 

6.63E-03 

3.25 

0.008 

2.27E-03 

4.08 

0.004 

2.49E-03 

4.07 

0.004 

40 

40 

2.54E-03 

1.39 

0.012 

l.OOE-04 

4.50 

0.016 

1.24E-04 

4.32 

0.016 

80 

80 

1.26E-03 

1.01 

0.028 

4.11E-06 

4.61 

0.036 

4.58E-06 

4.76 

0.040 

160 

160 

6.28E-04 

1.00 

0.124 

7.85E-07 

2.39 

0.124 

1.13E-07 

5.34 

0.152 

320 

320 

3.14E-04 

1.00 

0.424 

1.94E-07 

2.01 

0.464 

1.17E-08 

3.27 

0.528 

640 

640 

1.57E-04 

1.00 

1.668 

4.84E-08 

2.00 

1.805 

1.23E-09 

3.25 

2.128 


Table 6. Example 5 (ID Black and Scholes PDE). Error table 
with At ~ Ax, using SLDG-RKl, SLDG-RK2 and SLDG-RK3 
methods with P 4 polynomials (fc = 4). 


Example 6 (ID diffusion with non-constant cr(x)). Now, we consider the 
following diffusion equation 

Vt - ^cr'^{x)vj;a: = f{t, x), X G (0, 1), t € (0, T) 
x(0, x) = 0 X G (0,1), 
with periodic boundary conditions, 

(t(x) := sin( 27 rx). 


(113) 

(114) 
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and, for testing purposes, f{t,x) := Vt{t,x) — ^a‘^{x)vxx{t, x) where v{t,x) := 
sin(27rt) cos(27r(a; — t)), which is therefore the exact solution (v = v). 

In this case, in order to get higher than first-order accuracy in time, we use the 
SLDG-2 scheme corresponding to a Platen’s weak Taylor scheme. The correction 
for the source term f{t,x) is treated by adding the term (62) at Gauss quadrature 
points, at each time step. 

In Table we first check the accuracy with respect to time discretization, with 
fixed spatial mesh size so that only the time discretization error appears. 

Then, in Table the errors are given for varying mesh sizes such that At = Ax 
and with Pi or P 2 elements (fc = I or fc = 2). We find the expected orders for the 
schemes SLDG-1/2. 


Remark 5.2. Notice that there is no need for an assumption that the diffusion 
coefficient is non-vanishing in the proposed method. 


error 

SLDG-1 

SLDG-2 

N 

error 

order 

error 

order 

100 

1.19E-03 

- 

1.89E-04 

2.05 

200 

5.95E-04 

1.01 

4.57E-05 

1.97 

400 

2.96E-04 

1.01 

1.16E-05 

1.93 

800 

1.48E-04 

1.00 

3.07E-06 

1.91 

1600 

7.40E-05 

1.00 

8.17E-07 

1.92 


Table 7. Example 6 (ID diffusion with non-constant coefficient), 
with fixed spatial mesh {M = 100 and P 4 polynomials) and varying 
time steps N. 



error 

SLDG-1 (with Pi) 

SLDG-2 (with P 2 ) 

M 

N 

error 

order 

error 

order 

10 

10 

8.60E-02 

- 

4.13E-02 

- 

20 

20 

3.52E-02 

1.29 

7.30E-03 

2.50 

40 

40 

1.59E-02 

1.15 

1.39E-03 

2.39 

80 

80 

7.54E-03 

1.08 

3.03E-04 

2.20 

160 

160 

3.67E-03 

1.04 

7.17E-05 

2.08 

320 

320 

1.81E-03 

1.02 

1.80E-05 

1.99 


Table 8. Example 6 (ID diffusion with non-constant coefficient), 
with At ~ Ax; 


Example 7 (2D diffusion) We consider the following two-dimensional diffusion 
equation: 

Ut ^(.HUxx ^Uxy '^yy) — 0; X € D, t ^ (0, T), 
m(0,x) = Uq{x), X G 


(115) 

(116) 
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set on = (0,1)^ with periodic boundary conditions, and T = 0.2. The initial 
datum is given by uq{x) = Um{x + 2 y)+UQ 2 {-y) and uoi (0 := Y.q=i ,2 cos{2T:qC) 
with the constant c* = . The exact solution is knownj^ 

In order to define the numerical scheme, we use the fact that 

A := _2 j = X! 0-1 := f p j , (72 := 

The results are given in Table where we consider variable time steps and mesh 
steps At ^ Aa;, p = k, and expect a global error of order 0{At^) + ) = 

0{Ax'^). 

In this example involving constant diffusion coefficients, we test up to third-order 
schemes. 



error 

SLDG-RKl (Qi) 

SLDG-RK2 (Q 2 ) 

SLDG-RK3 (Q 3 ) 

II 

to 

N 

error 

order 

error 

order 

error 

order 

10 

10 

6.66E-03 

- 

1.86E-04 

- 

2 . 20 E -06 

- 

20 

20 

3.26E-03 

1.02 

4.52E-05 

2.04 

3.10E-07 

2.83 

40 

40 

1.61E-03 

1.01 

1.08E-05 

2.06 

3.20E-08 

3.27 

80 

80 

8.04E-04 

1.00 

2.69E-06 

2.01 

4.34E-09 

2.88 

160 

160 

4.01E-04 

1.00 

6.66E-07 

2.01 

4.90E-10 

3.14 


Table 9. Example 7 (2D diffusion equation), error table with 
At ~ Ax using Qu polynomials. 


Example 8 (2D diffusion with non-constant coefficients) We consider the 
following two-dimensional diffusion equation: 


ut - ]^Tr{aa'^D‘^u) = f{t,x), 


X e ft, t e (0,T), 


w(0, X, y) = uq{x, y), (x, y) 

set on n = (— 7 r, 7 r)^ with periodic boundary conditions, T = 1.0. 
matrix A = aa^ is defined by 


(117) 

(118) 
The diffusion 


a{x,y) := 


cos(x) cos(2x) 

0 sin(y) 

In this test we have chosen u{t, x, y) := cos(t)sin{2x)sin{x + y) and the source term 
f{t,x) such that (117) holds. (The initial datum is therefore uo{x,y) = u{0,x,y)). 

The scheme is defined here by using either 

• the weak Euler scheme for the diffusion part, combined with Trotter’s split¬ 
ting (and with Qi polynomials) and a first-order correction for the source 
tern (as in @)- 

• the weak Platen scheme for the diffusion part, combined with Strang’s 
splitting (and with Q 2 polynomials) and a second-order correction for the 
source term (roll) 


as explained in Section 4.2 


^ Making the change of variable J = (5i,?2) such that = x + 2y and ^2 = —y we find 
that v(t,() = u(t,x) satisfies vt — = 0 and = rtoi(Cl) + “ 02 ( 52 ) and 

therefore the exact solution is given by u(t,x) = u(t,5) = “i(i,5l) + “ 2 (i, 52 ) where Ui(b5) = 
^5=1,2 c* cos(27r55). 
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The results for errors are given in Table 10 where we consider variable time 
steps and mesh steps At ^ At. (see Section 4.1). The schemes are numerically 
roughly of the expected orders 1 and 2. 


As mentioned in Remark 1 5.2 [ there is no need to assume strict positivity of the 
diffusion matrix in this approach. 


error 

Euler/Trotter (with Qi) 

Platen/Strang (with Q 2 ) 

II 

A 

error 

order 

cpu(s) 

error 

order 

cpu(s) 

5 

10 

1.50E+00 

- 

0.01 

2.96E-01 

- 

0.020 

10 

20 

4.98E-01 

1.59 

0.02 

3.14E-02 

3.24 

0.088 

20 

40 

9.63E-02 

2.37 

0.11 

3.40E-03 

3.21 

0.432 

40 

80 

2.87E-02 

1.75 

0.74 

7.10E-04 

2.26 

2.564 

80 

160 

1.07E-02 

1.43 

5.44 

1.66E-04 

2.09 

16.621 


Table 10. Example 8 {2D diffusion equation with variable coefficients) 


Appendix A. Instability of the direct scheme 

Here we consider the ” direct scheme”, which defines naively at each time iteration 
a new piecewise polynomial € 14 such that, 

:= — 6At), for all Gauss points t),. 

In Figure we consider again Vt + = 0 with periodic boundary conditions 

on (0,1), and with the initial data vq{x) = sin(27rT). We have depicted two graphs 
with different choices of the parameter N. In each graph we plotted the result of 
the direct scheme (green line) and of the SLDG scheme (red line) at time T = 1, 
with piecewise Pi elements (fc = 1) and fixed spatial mesh using M = 46 mesh 
steps. In the left graph, A = 80 time steps and both curves are confounded; in the 
right graph, N = 320, and the direct scheme becomes unstable. (We have found 
that the error behaves as where c > 1, when using elements.) 




Figure 2. Results for A = 80 (left) and A = 320 (right), using 
Pi elements with M = 46 in both cases. Instability appears on the 
right figure. 
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